clear
rng default
addpath /Users/cristinagualdani/mosek/9.2/toolbox/r2015aom
%% Load useful matrices
load PXY_cond
PXY_cond_split=[PXY_cond(1:3); PXY_cond(4:end)];
load V_CS

%% Construct grid of parameters to check for each realisation of Y
gran=0.01;
v0grid_temp=0; %[Xchosen=0, Y=1] [Xchosen=0, Y=2], location normalisation 
v1grid_temp_all=[V_11_CS V_12_CS]; %[Xchosen=1, Y=1] [Xchosen=1, Y=2], scale normalisation
v2grid_temp=-20:gran:20;  %[Xchosen=2, Y=1] [Xchosen=2, Y=2]
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)

%% Construct identified set for each realisation of Y
ntypes_w=2;
IdSetW=cell(ntypes_w,1);


for y=1:ntypes_w %we can proceed separately across types as long as we do not impose independence from Y
    v1grid_temp=v1grid_temp_all(y);
    [ca, cb, cc] = ndgrid(v0grid_temp, v1grid_temp, v2grid_temp);
    v0grid=ca(:);
    v1grid=cb(:);
    v2grid=cc(:);
    sg=size(v0grid,1);
    % Components invariant across parameter values
    [equalorder, ...
    Ugridreduced, sgreduced,...
    id_Ugridreduced,...
    s_id_gridreduced, sU_id_gridreduced,...
    n_U_sets, Tcoord, indices_pairs,...
    numeqconstraints,coordrowseq, fillAeq]=useful_anyparam4a(v0grid, v1grid, v2grid, sg);
    % Loop over parameter values, order of columns: [u0x u1x u2x]
    IdSetW{y}=OneSide4a(PXY_cond_split(y,:), v0grid, v1grid, v2grid,...
                        equalorder, ...
                        Ugridreduced, sgreduced,...
                        id_Ugridreduced,...
                        s_id_gridreduced, sU_id_gridreduced,...
                        n_U_sets, Tcoord, indices_pairs,...
                        numeqconstraints,coordrowseq, fillAeq);  
end


[cx1, cx2] = ndgrid(1:size(IdSetW{1},1), 1:size(IdSetW{2},1));
cx1=cx1(:);
cx2=cx2(:);
IdSetW_final=[IdSetW{1}(cx1,:) IdSetW{2}(cx2,:)]; 
%order of columns: [V01 V11 V21 V02 V12 V22]

save('IdSetW_final.mat', 'IdSetW_final')